Planejamento experimental
e Análise de dados agronômicos

José Tiago Barroso Chagas

Apresentação

Engenheiro agronômo (UFCA- 2018)

Mestre em genética e melhoramento de plantas (UENF-2020)

Doutor em genética e melhoramento (UFV- 2024)

Pesquisador de pós doutorado (UFV- Atual)

Sumário

  • Princípios da experimentação agrícola
  • Planejamento de Experimentos agronômicos
  • Equações de Modelos Mistos
  • Análise de variância individual
  • Análise de variância conjunta
  • Predição de valores genotípicos
  • Análise de Fatores
  • Índices de seleção

Princípios da Experimentação agrícola

  • Repetição
  • Casualização
  • Controle local
DIC - 5trat*3rep

DBC - 5trat*3rep

Planejamento de Experimentos agronômicos

  • Experimentos com grande número de tratamentos para os blocos casualizados em DBC
  • k Tratamentos divididos em k blocos (ktrat*kblocos)
  • Látice quadrado (9trat/k3/rep3)
LÁTICE

  ID LOCATION PLOT REP IBLOCK UNIT ENTRY TREATMENT
1  1  COIMBRA 1001   1      1    1     7       G-7
2  2  COIMBRA 1002   1      1    2     3       G-3
3  3  COIMBRA 1003   1      1    3     4       G-4

Planejamento de Experimentos agronômicos

  • Experimento com restrição quanto à disponibilidade de sementes e/ou volume grande de tratamentos para a área experimental
  • Delineamento em blocos aumentados

(30gen/3test/k6)

  ID  EXPT LOCATION YEAR PLOT ROW COLUMN CHECKS BLOCK ENTRY TREATMENT
1  1 Expt1  COIMBRA 2024 1001   1      1      1     1     3       CH3
2  2 Expt1  COIMBRA 2024 1002   1      2      1     1     2       CH2
3  3 Expt1  COIMBRA 2024 1003   1      3      0     1    23       G23
DBA

::::

Equações de modelos mistos

Charles R. Henderson (1953) (melhorista animal) propôs as equações de modelos mistos

As quais possibilitaram obter BLUPs dos efeitos aleatórios (Melhor preditor linear não viesado) e estimadores de máxima verossimilhança dos efeitos fixos de \(\beta\).

Equações de modelos mistos

\[ \textbf{y} = \mathbf{X}\beta + \mathbf{Z}u+ \mathbf{e} \]

\[ \begin{align*} \begin{bmatrix} \mathbf{X}'\mathbf{X} & \mathbf{X}'\mathbf{Z} \\ \mathbf{Z}'\mathbf{X} & \mathbf{Z}'\mathbf{Z} + \lambda \mathbf{K}^{-1} \end{bmatrix} \begin{bmatrix} \mathbf{\hat{b}} \\ \mathbf{\hat{u}} \end{bmatrix} &= \begin{bmatrix} \mathbf{X}'\mathbf{y} \\ \mathbf{Z}'\mathbf{y} \end{bmatrix} \end{align*} \]

Máxima verossimilhança restrita
(REML-Restricted maximum likelihood)

Patterson and Thompson (1971)

Função de densidade de probabilidade normal

\[ f(x)= \frac{1}{\sigma \sqrt{2\pi}}e^{-\frac{1}{2}(\frac{x-\mu}{\sigma})^2} \]

A função de máxima verossimilhança buscar a otimização dos valores da média \(\mu\) e variância \(\sigma\), dado as observações e considerando distribuição normal.

Neste exemplo consideramos 100 valores aleatórios gerados de uma distribuição normal com média 10 (\(\mu\)) e desvio padrão \(\sqrt{10}\) (\(\sigma\))

head de valores amostrados
[1] 13.105260 11.482207  9.658565  9.326820 13.662229 14.086785

Equações de modelos mistos

BLUE (Melhor estimador linear não viesado - efeitos fixos)

BLUP (Melhor preditor linear não viesado - efeitos aleatórios)

C. R. Henderson (1975)

Simulação de dados fenotípicos - Exemplo

Traits: Produtividade de grãos e altura de planta

Fonte: Bayer crop science

Fonte: Bayer crop science

2 características; 100 genótipos; 2 blocos

Variâncias residuais Trait per environment -
    c(Yld:E1, Yld:E2, Yld:E3, Pht:E1, Pht:E2, Pht:E3
[1]  0.20  0.28  0.14 15.10  8.50 11.70
Correlação dos erros devido tendências espaciais
     1    2
1 1.00 0.49
2 0.49 1.00
Correlação dos erros devido causas não conhecidas
    (Caminhamento,máquinas e implementos)
     1    2
1 1.00 0.02
2 0.02 1.00
Banco de dados simulado de Milho
'data.frame':   200 obs. of  8 variables:
 $ env  : Factor w/ 1 level "1": 1 1 1 1 1 1 1 1 1 1 ...
 $ block: Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
 $ col  : Factor w/ 10 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ row  : Factor w/ 20 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
 $ id   : Factor w/ 100 levels "1","2","3","4",..: 39 4 77 14 49 62 47 94 40 63 ...
 $ rep  : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
 $ GY   : num  5.51 5.7 5.8 5.17 5.32 ...
 $ PH   : num  230 230 230 232 239 ...

Solving mixed models with RHS
[1] 5.513532 5.703086 5.796593 5.168647 5.323202 5.033188

Resolvendo as equações de modelos mistos (RHS) - GY

\[ \begin{align*} \begin{bmatrix} \mathbf{X}'\mathbf{X} & \mathbf{X}'\mathbf{Z} \\ \mathbf{Z}'\mathbf{X} & \mathbf{Z}'\mathbf{Z} + \lambda \mathbf{K}^{-1} \end{bmatrix} \begin{bmatrix} \mathbf{\hat{b}} \\ \mathbf{\hat{u}} \end{bmatrix} &= \begin{bmatrix} \mathbf{X}'\mathbf{y} \\ \mathbf{Z}'\mathbf{y} \end{bmatrix} \end{align*} \]

X matrix - Efeitos fixos
  block1 block2
1      1      0
2      1      0
3      1      0
4      1      0
5      1      0
6      1      0
Z matrix - Efeitos aleatórios
   id1 id2 id3 id4 id5 id6 id7 id8 id9 id10
1    0   0   0   0   0   0   0   0   0    0
2    0   0   0   1   0   0   0   0   0    0
3    0   0   0   0   0   0   0   0   0    0
4    0   0   0   0   0   0   0   0   0    0
5    0   0   0   0   0   0   0   0   0    0
6    0   0   0   0   0   0   0   0   0    0
7    0   0   0   0   0   0   0   0   0    0
8    0   0   0   0   0   0   0   0   0    0
9    0   0   0   0   0   0   0   0   0    0
10   0   0   0   0   0   0   0   0   0    0
Phenotypic data
   env block col row id rep       GY       PH
1    1     1   1   1 39   1 5.513532 230.2357
2    1     1   1   2  4   1 5.703086 230.2436
3    1     1   1   3 77   1 5.796593 230.4429
4    1     1   1   4 14   1 5.168647 231.6617
5    1     1   1   5 49   1 5.323202 238.9142
6    1     1   1   6 62   1 5.033188 233.4361
7    1     1   1   7 47   1 5.969180 227.9464
8    1     1   1   8 94   1 5.296897 233.5849
9    1     1   1   9 40   1 3.751499 226.4649
10   1     1   1  10 63   1 4.588704 235.6345
Error variance
[1] 0.2075
Genotypic variance
[1] 0.0948
Lambda
[1] 2.188819
XpX
       block1 block2
block1    100      0
block2      0    100
XpZ
       id1 id2 id3 id4 id5
block1   1   1   1   1   1
block2   1   1   1   1   1
ZpX
    block1 block2
id1      1      1
id2      1      1
id3      1      1
id4      1      1
id5      1      1
id6      1      1
ZpZ
    id1 id2 id3 id4 id5
id1   2   0   0   0   0
id2   0   2   0   0   0
id3   0   0   2   0   0
id4   0   0   0   2   0
id5   0   0   0   0   2
Xpy
           [,1]
block1 498.0205
block2 509.2963
Zpy
         [,1]
id1  8.834306
id2 10.478144
id3 11.530348
id4 10.042302
id5  9.734658
id6  9.007352
LHS
       block1 block2      id1      id2      id3
block1    100      0 1.000000 1.000000 1.000000
block2      0    100 1.000000 1.000000 1.000000
id1         1      1 4.188819 0.000000 0.000000
id2         1      1 0.000000 4.188819 0.000000
id3         1      1 0.000000 0.000000 4.188819
RHS
             [,1]
block1 498.020468
block2 509.296324
id1      8.834306
id2     10.478144
id3     11.530348
id4     10.042302
Solutions solve(LHS) %*% RHS
               [,1]
block1  4.980204681
block2  5.092963242
id1    -0.295754501
id2     0.096680330
id3     0.347873860
id4    -0.007368719

Resolvendo as equações de modelos mistos (LME4) - GY

Solving mixed models with package LME4
Linear mixed model fit by REML ['lmerMod']
Formula: GY ~ rep + (1 | id)
   Data: pheno_env1

REML criterion at convergence: 320.2

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.15158 -0.56334  0.05882  0.58978  2.09491 

Random effects:
 Groups   Name        Variance Std.Dev.
 id       (Intercept) 0.1010   0.3178  
 Residual             0.1982   0.4451  
Number of obs: 200, groups:  id, 100

Fixed effects:
            Estimate Std. Error t value
(Intercept)  4.98020    0.05469  91.058
rep2         0.11276    0.06295   1.791

Correlation of Fixed Effects:
     (Intr)
rep2 -0.576
RHS
x
-0.2957545
0.0966803
0.3478739
-0.0073687
-0.0808127
-0.2544430
LME4
x
id.(Intercept)1 -0.3126610
id.(Intercept)2 0.1022069
id.(Intercept)3 0.3677597
id.(Intercept)4 -0.0077899
id.(Intercept)5 -0.0854323
id.(Intercept)6 -0.2689879

Análise conjunta

\[ \mathbf{y} = \mathbf{X}_1\mathbf{e} + \mathbf{X}_2\mathbf{be} + \mathbf{Z}_1\mathbf{g} + \mathbf{Z}_2\mathbf{ge} + \mathbf{\epsilon} \]

\([\mathbf{g}\sim N(0,\sigma_g^2\mathbf{I}_n)]\) e \([\mathbf{ge}\sim N(0,\sigma_{ge}^2\mathbf{I}_{np})]\), onde n é o número de genótipos e p é o número de ambientes.

2 características; 100 genótipos; 2-3 blocos; 3 ambientes

Testar os efeitos de genótipo e da interação genótipo x ambiente para as características avaliadas. \[ LRT = -2(LogL_-LogL_{r}) \] onde \(LogL\) é o ponto de máxima da função de verossimilhança residual do modelo completo e \(LogL_{r}\) é o mesmo para o modelo reduzido.

Análise de variância conjunta - GY
Linear mixed model fit by REML ['lmerMod']
Formula: GY ~ env + rep:env + (1 | id) + (1 | id:env)
   Data: pheno_df

REML criterion at convergence: 1063.2

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.67400 -0.55510  0.01111  0.57905  2.30558 

Random effects:
 Groups   Name        Variance Std.Dev.
 id:env   (Intercept) 0.09360  0.3059  
 id       (Intercept) 0.01043  0.1021  
 Residual             0.18065  0.4250  
Number of obs: 700, groups:  id:env, 300; id, 100

Fixed effects:
            Estimate Std. Error t value
(Intercept)  4.98020    0.05336  93.340
env2         0.13322    0.07406   1.799
env3         0.48945    0.07406   6.609
env1:rep2    0.11276    0.06011   1.876
env2:rep2   -0.02553    0.06011  -0.425
env3:rep2   -0.12331    0.06011  -2.051
env3:rep3   -0.19016    0.06011  -3.164

Correlation of Fixed Effects:
          (Intr) env2   env3   env1:2 env2:2 env3:2
env2      -0.694                                   
env3      -0.694  0.500                            
env1:rep2 -0.563  0.406  0.406                     
env2:rep2  0.000 -0.406  0.000  0.000              
env3:rep2  0.000  0.000 -0.406  0.000  0.000       
env3:rep3  0.000  0.000 -0.406  0.000  0.000  0.500
fit warnings:
fixed-effect model matrix is rank deficient so dropping 2 columns / coefficients
LRT GxE effect - Trait GY
[1] 43.31988
Data: pheno_df
Models:
mmcr: GY ~ env + rep:env + (1 | id)
mmc: GY ~ env + rep:env + (1 | id) + (1 | id:env)
     npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
mmcr    9 1096.1 1137.0 -539.03   1078.1                         
mmc    10 1054.3 1099.8 -517.15   1034.3 43.758  1  3.717e-11 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
LME4 - g x e
(Intercept)
id:env.1:1 -0.2994339
id:env.1:2 -0.0883629
id:env.1:3 0.1092676
id:env.2:1 0.0992936
id:env.2:2 -0.0491773
id:env.2:3 0.0160696
LME4 - g
(Intercept)
id.1 -0.0310489
id.2 0.0073780
id.3 0.0317793
id.4 -0.0026264
id.5 0.0316266
id.6 -0.0007239

Índice de Seleção

Diversos índices de seleção foram propostos com e sem restrições (Smith (1936) & Hazel (1943) ; Kempthorne and Nordskog (1959))

Neste curso abordaremos o índice de seleção da análise de fatores e desenho da distância genótipo ideótipo Rocha, Machado, and Carneiro (2018).

Análise de Fatores

Considerando o dataset simulado de milho (2 traits * 3 ambientes)

Para proceder análise fatorial, deve-se realizar a análise de correlação de pearson entre pares de variáveis e criar a matriz de correlações.

Blups per trait:env
          GY_1        GY_2         GY_3       PH_1       PH_2       PH_3
1 -0.299433931 -0.08836286  0.109267596 -0.7960276 -3.2359131  1.5381855
2  0.099293639 -0.04917734  0.016069572 -2.6922993 -2.2370177  1.3938318
3  0.354614536  0.04186552 -0.111398567  0.7017383 -0.4216252  1.2066291
4 -0.006517491 -0.23820750  0.221164735 -3.2523328 -1.3943519 -0.1705401
5 -0.102230756  0.26035403  0.125588089 -1.4563768 -3.3728975  1.4816980
6 -0.270834073  0.26110773  0.003232461 -7.0252473 -0.8419086  2.4660021
Blups per trait:env standarized
         GY_1       GY_2        GY_3       PH_1       PH_2        PH_3
1 -1.37867867 -0.3996783  0.49698215 -0.4213632 -1.5317044  0.80693911
2  0.45717605 -0.2224364  0.07308928 -1.4251213 -1.0588819  0.73121050
3  1.63274581  0.1893639 -0.50667445  0.3714528 -0.1995743  0.63300313
4 -0.03000838 -1.0774477  1.00592425 -1.7215651 -0.6600100 -0.08946613
5 -0.47069937  1.1776197  0.57121270 -0.7709074 -1.5965453  0.77730550
6 -1.24699682  1.1810288  0.01470221 -3.7186910 -0.3985135  1.29367592

Matriz de correlações

            GY_1        GY_2        GY_3        PH_1       PH_2        PH_3
GY_1  1.00000000 -0.05369986 -0.15382851  0.09624741  0.1535808 -0.03312696
GY_2 -0.05369986  1.00000000  0.04461362 -0.08982962  0.2295911  0.14690621
GY_3 -0.15382851  0.04461362  1.00000000  0.09267733 -0.2505011  0.38971600
PH_1  0.09624741 -0.08982962  0.09267733  1.00000000  0.1306623 -0.26498653
PH_2  0.15358083  0.22959111 -0.25050112  0.13066234  1.0000000 -0.34163029
PH_3 -0.03312696  0.14690621  0.38971600 -0.26498653 -0.3416303  1.00000000

Correlações próximas de zero propiciam a extração de fatores diferentes. Representar as correlações entre os traits por meio de variáveis latentes (fatores).

Autovalores

Determinação dos autovalores pela seguinte equação:

\[ [\det(\lambda^².I-\rho) = 0] \]

Multiplicado por uma matriz identidade diminuida da matriz de correlações \(\rho\)

\[ \begin{bmatrix} \lambda^2 & 0 & 0 & \cdots & 0 \\ 0 & \lambda^2 & 0 & \cdots & 0 \\ 0 & 0 & \lambda^2 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & \lambda^2 \end{bmatrix} \begin{bmatrix} \rho_{11} & \rho_{12} & \rho_{13} & \cdots & \rho_{1n} \\ \rho_{21} & \rho_{22} & \rho_{23} & \cdots & \rho_{2n} \\ \rho_{31} & \rho_{32} & \rho_{33} & \cdots & \rho_{3n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \rho_{n1} & \rho_{n2} & \rho_{n3} & \cdots & \rho_{nn} \end{bmatrix}\cdots \begin{bmatrix} \lambda^2 - \rho_{11} & -\rho_{12} & -\rho_{13} & \cdots & -\rho_{1n} \\ -\rho_{21} & \lambda^2 - \rho_{22} & -\rho_{23} & \cdots & -\rho_{2n} \\ -\rho_{31} & -\rho_{32} & \lambda^2 - \rho_{33} & \cdots & -\rho_{3n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ -\rho_{n1} & -\rho_{n2} & -\rho_{n3} & \cdots & \lambda^2 - \rho_{nn} \end{bmatrix}= 0 \]

Polinômio característico

Obtenção do polinôminio característica de ordem relativa à quantidade de fatores.

(n características - ordem n)

Error term: 1.070302 
Polinômio característico
0.5352 - 3.9883*x + 11.8003*x^2 - 17.7752*x^3 + 14.4279*x^4 - 6*x^5 + x^6 
Raízes do polinômio - Autovalores
[1] 1.7653 1.2008 1.0624 0.9772 0.5603 0.4340

Os auto valores são as raízes do polinômio característico da matriz de correlações.

A soma dos autovalores é igual a quantidade de variáveis originais

[1] 6

Percentual de variância - valor do autovalor dividido pelo total

[1] 29.4224 20.0130 17.7060 16.2864  9.3389  7.2333

Autovetores

Os autovetores são derivados dos autovalores e resultam de um sistema de equações.

\[ [det(\lambda^².I-\rho)\mathbf{v}_i = 0] \]

\[ \begin{bmatrix} \lambda^2 - \rho_{11} & -\rho_{12} & -\rho_{13} & \cdots & -\rho_{1n} \\ -\rho_{21} & \lambda^2 - \rho_{22} & -\rho_{23} & \cdots & -\rho_{2n} \\ -\rho_{31} & -\rho_{32} & \lambda^2 - \rho_{33} & \cdots & -\rho_{3n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ -\rho_{n1} & -\rho_{n2} & -\rho_{n3} & \cdots & \lambda^2 - \rho_{nn} \end{bmatrix} \begin{bmatrix} v_{i1} \\ v_{i2} \\ v_{i3} \\ \vdots \\ v_{in} \end{bmatrix}= 0 \]

\[ (\lambda^2 - \rho_{11}) v_{i1} - \rho_{12} v_{i2} - \rho_{13} v_{i3} - \cdots - \rho_{1n} v_{in} = 0 ] \]

Escores fatoriais

Obtido por meio do autovalor dividido pela raiz quadrada do seu respectivo autovetor.

Indica relação com as variáveis originais.

\[ S_{11}=\frac{v_{11}}{\sqrt{\lambda²_1}} \]

Extração dos fatores

Os escores fatoriais multiplicados pelas variáveis padronizadas (Zscore) geram os fatores.

A grandeza dos escores para cada variável representará o seu peso em cada fator.

\[ S_{11}=\frac{v_{11}}{\sqrt{\lambda²_1}} \]

Cargas fatoriais

Representam as correlações de pearson entre os fatores extraídos e as variáveis originais (+carga = +importância da variável em cada fator).

O resultado da soma das cargas fatoriais ao quadrado de cada fator é próprio autovalor do qual ele foi extraído.

Mais informações sobre a história da análise fatorial nesse artigo Bartholomew (1995)

Índice de Seleção

\[ Pij = \frac{\frac{1}{d_{ij}}}{\sum^{i=n;j=m}_{i=1;j=1} \frac{1}{d_{ij}}} \]

\(Pij\) probabilidade do i-ésimo genótipo ser similar ao j-ésimo ideótipo, \(d_{ij}\) distância genótipo-ideótipo euclidiana padronizada.

Referências

Bartholomew, D. J. 1995. “Spearman and the Origin and Development of Factor Analysis.” British Journal of Mathematical and Statistical Psychology 48 (November): 211–20. https://doi.org/10.1111/J.2044-8317.1995.TB01060.X.
Hazel, L N. 1943. “The Genetic Basis for Constructing Selection Indexes.” Genetics 28: 476–90. http://www.ncbi.nlm.nih.gov/pubmed/17247099%0Ahttp://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=PMC1209225.
Henderson, C R. 1975. “Best Linear Unbiased Estimation and Prediction Under a Selection Model.” Biometrics 31: 423. https://doi.org/10.2307/2529430.
Henderson, Charles R. 1953. “Estimation of Variance and Covariance Components.” Biometrics 9 (2): 226–52.
Kempthorne, Oscar, and Arne W Nordskog. 1959. “Restricted Selection Indices.” Biometrics 15: 10–19. http://www.jstor.org/stable/2527598.
Patterson, H D, and R Thompson. 1971. “Recovery of Inter-Block Information When Block Sizes Are Unequal.” Biometrika 58: 545–54. https://doi.org/10.2307/2334389.
Rocha, João Romero do Amaral Santos de Carvalho, Juarez Campolina Machado, and Pedro Crescêncio Souza Carneiro. 2018. “Multitrait Index Based on Factor Analysis and Ideotype-Design: Proposal and Application on Elephant Grass Breeding for Bioenergy.” Gcb Bioenergy 10 (1): 52–60.
Smith, H. F. 1936. “A Discriminant Function for Plant Selection.” Annuals of Eugenics 7: 240–50.

Exp_Plan_&_Data_Analysis